% rsts_STAT.M
% compute posterior mean, sd, and confidence interval (HPD)
% based on Metropolis-Hastings output

% Taeyoung Doh
% created: 02/02/2007
%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
loaddata;

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);


nobs    = size(data,1); 

global datasel psel 
% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for linear posterior draws
%	even integer (e.g., 2) for quadratic posterior draws
 
mhrun = 1;
mhrunid  = num2str(mhrun','%2.0f');

block1   = 1;		% starting block
blockn   = 150;		% ending block
nsim     = 10000;	% number of simulation draws in each block

%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% parameter names
pnames = feval(fnames_.pnames,msel);

% reset diary file
diaryfile = [path_.result 'stat_mhrun' mhrunid '.log'];
if exist(diaryfile,'file')
	delete(diaryfile);
end
diary(diaryfile);
diary off;

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


% load path
  lpath = [path_.para 'mhrunpmchain2' mhrunid '/' ];  



% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);
    [B,I]=sort(postsim); 
    if jblock==block1
    drawmax=parasim(I(end),:); 
    postmax=B(end);
    likemax=likesim(I(end));
    else 
           if postmax<B(end)
              drawmax=parasim(I(end),:);
              postmax=B(end);
              likemax=likesim(I(end));
           end
    end 
    

    
    
	
end
drawdim=length(drawmax);
%--------------------------------------------------------------------
% save and print
%--------------------------------------------------------------------
% save output
%ssfile = [path_.para 'rs_prob' mhrunid '.mat'];
%save(ssfile,'ssT'); 

% print log
fmt='%15.8f ';
disp(' ');
disp(' ');
disp(' ');
diary on;
disp(' ');
disp(sprintf('-------------------------------------------------------------------------'));
disp(         'Parameters     Estimates ');
for i=1:drawdim 
	disp(sprintf(['%s '  repmat(fmt,1,2) ],pnames(i,:),drawmax(1,i)));
end
disp(sprintf('-------------------------------------------------------------------------'));
diary off;


sfile = [path_.para 'pmode.mat'];
save(sfile,'drawmax','postmax','likemax');

